Visualisations of a species’ distribution can be some of the most simple yet effective conveyors of biological and ecological information, such as range, habitat, population size, fragmentation and migrations. Key to their efficacy is their versatility and variability - distributions can be clearly represented as points, hulls, densities, and binned regions, to name a few.
However, visualising more than one species’ distribution is more difficult to perform on a single map, especially when those distributions overlap. Points and hulls will obscure each other even with a degree of transparency, while densities and shaded regions can only show one species at a time. Comparisons can be made by providing two single distribution figures side-by-side, but it can be inconvenient to observe range overlaps in this manner.
Here, we will demonstrate a method to visualise distributions of multiple species with overlapping range on the same map, for a small loss in resolution. The technique is a novel twist on the commonly used hex maps, - instead of shading we use coloured (or shaped) points within the hexagons to indicate the presence or absence of a species.
Download Data
We begin by loading necessary R packages.
library(galah)
library(ggtext)
library(glue)
library(hexbin)
library(monochromeR)
library(ozmaps)
library(sf)
library(tidyverse)We will use the {galah} package to download occurrence records from the Atlas of Living Australia (ALA) for some chosen group of species. Make sure you provide a registered email address with the ALA through galah_config(). The default settings provided by galah_config() are sufficient for what we wish to do.
galah_config(email = "your-email@email.com")Download Shapefiles
If we are going to visualise species’ distributions then we need a relevant base map to plot those distributions on top of. In this case, we are going to begin with all land that falls under Australia’s jurisdiction (including islands such as Christmas Island, Norfolk Island, Macquarie Island etc.). The base layer we will use for plotting will come from the {ozmaps} package which provides simple shapefiles of terrestrial Australia. We can save this object as aus.
# map of Australia
aus <- ozmap_data(data = "states")However, higher resolution shapefiles exist of Australian land, and one such resource is the Interim Biogeographic Regionalisation of Australia (IBRA). The shapefile and related metadata can be downloaded from this DCCEEW link. Because of it’s finer detail, the IBRA dataset is better suited for filtering records that occurred within Australia than the {ozmaps} layer, but does not contain state boundaries which are useful for plotting. We will use both maps of Australia throughout our workflow and visualisation.
Once you’ve downloaded and saved this shapefile in your working directory, we use the {sf} package to read the file into R. Any functions we use in the workflow that begin with the prefix st_ are from {sf} and are used primarily for manipulation of spatial objects. We call st_make_valid() within the pipe to ensure that we have a set of valid polygons, and then st_union() to create a single spatial feature from all of the IBRA regions that does not contain any internal borders. This code block may take a little bit of time to run so be patient.
ibra <- st_read("path/to/IBRA7_regions/ibra7_regions.shp") |>
st_make_valid() |>
st_union()
crs <- st_crs(ibra)We also store the coordinate reference system (CRS) of the IBRA layer using st_crs(), which we will use later on to ensure all our sf objects follow the same CRS.
Download Occurrence Data
We next must choose a group of species to visualise. Our goal here is to map distributions of occurrences of more than one species, and here we are going to do it for 7. In this case we have chosen the black-headed honeyeaters, genus Melithreptus, which are abundant across Australia. We can verify that there are indeed 7 such species in the ALA using atlas_species().
galah_call() |>
galah_identify("Melithreptus") |>
atlas_species()# A tibble: 7 × 10
kingdom phylum class order family genus species author species_guid
<chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… (Viei… https://bio…
2 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… (Vigo… https://bio…
3 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… Gould… https://bio…
4 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… (Goul… https://bio…
5 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… (Less… https://bio…
6 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… (Goul… https://bio…
7 Animalia Chordata Aves Passeriformes Melip… Meli… Melith… Gould… https://bio…
# ℹ 1 more variable: vernacular_name <chr>We can create a simple dataframe, species_data as follows that contains the species names, common names, and a colour that we will assign to each species.
species_data <- data.frame(
species = c("Melithreptus albogularis", # White-throated Honeyeater
"Melithreptus lunatus", # White-naped Honeyeater
"Melithreptus brevirostris", # Brown-headed Honeyeater
"Melithreptus gularis", # Black-chinned Honeyeater
"Melithreptus affinis", # Black-headed Honeyeater
"Melithreptus chloropsis", # Gilbert's Honeyeater
"Melithreptus validirostris"), # Strong-billed Honeyeater
common_name = c("White-throated Honeyeater",
"White-naped Honeyeater",
"Brown-headed Honeyeater",
"Black-chinned Honeyeater",
"Black-headed Honeyeater",
"Gilbert's Honeyeater",
"Strong-billed Honeyeater"),
colour = c("#F7618F",
"#842192",
"#F7C328",
"#33C8E1",
"#E4C9C9",
"#D7271C",
"#7CC545")
) |>
mutate(label = paste(common_name, " (*", species, "*)", sep = ""))
species_data species common_name colour
1 Melithreptus albogularis White-throated Honeyeater #F7618F
2 Melithreptus lunatus White-naped Honeyeater #842192
3 Melithreptus brevirostris Brown-headed Honeyeater #F7C328
4 Melithreptus gularis Black-chinned Honeyeater #33C8E1
5 Melithreptus affinis Black-headed Honeyeater #E4C9C9
6 Melithreptus chloropsis Gilbert's Honeyeater #D7271C
7 Melithreptus validirostris Strong-billed Honeyeater #7CC545
label
1 White-throated Honeyeater (*Melithreptus albogularis*)
2 White-naped Honeyeater (*Melithreptus lunatus*)
3 Brown-headed Honeyeater (*Melithreptus brevirostris*)
4 Black-chinned Honeyeater (*Melithreptus gularis*)
5 Black-headed Honeyeater (*Melithreptus affinis*)
6 Gilbert's Honeyeater (*Melithreptus chloropsis*)
7 Strong-billed Honeyeater (*Melithreptus validirostris*)# 6 5
# 1 7 4
# 2 3With a set of species names, we can now download occurrences of these species. As this is a bird genus, it should be no surprise that there are hundreds of thousands of records in the ALA (>516,000). These would take a very long time to download so here we will just download occurrences from 2022, of which there are ~35,500.
Using {galah}, we can obtain all such occurrences and store them as species_occ. We filter to obtain only records from 2022 for taxa in the species_data$species column. We also provide additional filters on occurrenceStatus, location and date parameters to ensure we are only keeping complete occurrences of birds. The final filter applies to the IBRA regions contextual layer, cl1048, and removes any records that did not occur in an IBRA region.
For each occurrence, we also request additional information on top of the basis info, including species, vernacularName, and cl1048 (IBRA Regions layer). We specifically request the species field with galah_select(), as the scientificName field often contains subgenera and subspecies names. We do not want this much detail - we simply want to match species names to our species_data$species column.
species_occ <- galah_call() |>
galah_apply_profile(ALA) |>
galah_filter(year == 2022) |>
galah_filter(occurrenceStatus == "PRESENT",
decimalLatitude != "",
decimalLongitude != "",
eventDate != "",
cl1048 != "") |>
galah_identify(species_data$species) |>
galah_select(group = c("basic"), species, vernacularName, cl1048) |>
atlas_occurrences()Request for 35550 occurrences placed in queue
Current queue length: 1---Downloadinghead(species_occ)# A tibble: 6 × 11
recordID scientificName taxonConceptID decimalLatitude decimalLongitude
<chr> <chr> <chr> <dbl> <dbl>
1 000096cc-6ec8-… Melithreptus … https://biodi… -27.3 153.
2 0001625a-18f5-… Melithreptus … https://biodi… -13.0 131.
3 0003ba50-1c17-… Melithreptus … https://biodi… -43.0 147.
4 00055052-798b-… Melithreptus … https://biodi… -35.5 149.
5 0006f874-c8ce-… Melithreptus … https://biodi… -16.6 145.
6 0009d7a1-8dd2-… Melithreptus … https://biodi… -34.2 116.
# ℹ 6 more variables: eventDate <dttm>, occurrenceStatus <chr>,
# dataResourceName <chr>, species <chr>, vernacularName <chr>, cl1048 <chr>Introduce Spatial Features
Occurrence Locations
Our species_occ dataframe contains all records from 2022 of the genus Melithreptus that we can use to identify species distributions. The next step is to convert our data frame into an sf object with a geometry column representing the location of each occurrence. We use st_as_sf() to convert the decimalLatitude and decimalLongitude variables into a single column of spatial POINT objects. We use the stored IBRA CRS (crs) and rename the geometry column as occ_sf using st_geometry(). In this workflow, column names ending in _sf will denote those containing spatial objects.
species_occ_sf <- species_occ |>
st_as_sf(coords = c("decimalLongitude", "decimalLatitude"),
crs = st_crs(crs), remove = FALSE)
st_geometry(species_occ_sf) <- "occ_sf"Hexagons
With the species data downloaded, the next step is to set up the hex grid over the IBRA regions. We can use st_make_grid() from {sf} to easily perform this. Note the variety of arguments in this function that provide flexibility - varying cellsize controls the “diameter” of the hexagons, flat_topped controls the orientation of the hexagons, and square enables the creation of a square grid instead, which may be more suitable for certain numbers of species (4,5,8,9 etc). We also assign a unique ID number to each hexagon, and rename the spatial column hex_sf.
all_hex_grid <- st_make_grid(ibra,
cellsize = 2,
what = "polygons",
square = FALSE,
flat_topped = TRUE) |>
st_as_sf() |>
mutate(hex_id = 1:n())
st_geometry(all_hex_grid) <- "hex_sf"As this creates a complete grid across the smallest rectangle containing all Australian IBRA regions, we next must restrict our hexagons to only those that actually contain Australian land. First, we use st_intersects() to extract the IDs of the hexagons that overlap with our ibra shapefile. We can then filter all_hex_grid to keep those hexagons. Because we’ll need it later, we’ll also take the opportunity in the pipe to add a column, hex_centre_sf, to our new dataframe ibra_hex_grid that contains the centre of each hexagon using st_centroid().
ibra_hex_ids <- st_intersects(all_hex_grid, ibra) |>
as.data.frame() |>
pull(row.id)
ibra_hex_grid <- all_hex_grid |>
filter(hex_id %in% (ibra_hex_ids)) |>
mutate(hex_centre_sf = st_centroid(hex_sf))
head(ibra_hex_grid)Simple feature collection with 6 features and 1 field
Active geometry column: hex_sf
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 71.42268 ymin: -55.77699 xmax: 160.3346 ymax: -41.77699
Geodetic CRS: GDA94
hex_sf hex_id hex_centre_sf
1 POLYGON ((158.0252 -54.7769... 54 POINT (159.1799 -54.77333)
2 POLYGON ((73.15473 -53.7769... 57 POINT (74.30943 -53.77355)
3 POLYGON ((71.42268 -52.7769... 85 POINT (72.57738 -52.77376)
4 POLYGON ((145.9009 -43.7769... 358 POINT (147.0556 -43.77521)
5 POLYGON ((144.1688 -42.7769... 386 POINT (145.3235 -42.77534)
6 POLYGON ((147.6329 -42.7769... 387 POINT (148.7876 -42.77534)Our map of hexagons now look like this,
ggplot() +
geom_sf(data = aus,
col = "grey20", fill = NA, linewidth = 0.3) +
geom_sf(data = ibra,
col = "grey20", fill = NA, linewidth = 0.3) +
geom_sf(data = ibra_hex_grid, aes(geometry = hex_sf),
fill = NA, col = "grey60", linewidth = 0.5) +
theme_void()
but it should be clear that some of these hexagons won’t contain any occurrences of Melithreptus. We will want to remove those hexagons that don’t contain any records to produce a cleaner map.
Species and Hexagon Combinations
To do this we have to find all unique combinations of species and hexagons by identifying which hexagon each occurrence is located in. To do this we use a neat application of mutate() and st_intersects() as described in this helpful GIS StackExchange answer. Taking our data-frame of species occurrences , we create a column intersection that returns the row index of ibra_hex_grid of the hexagon the record was located. In the same mutate() call, we then bring in the ID of that hexagon by calling the hex_id of that row index in intersection. We then only wish to retain unique combinations of species and hexagon IDs so we drop the spatial feature (st_drop_geometry() from the new data-frame and then perform the required operations
species_hex <- species_occ_sf |>
mutate(intersection = st_intersects(occ_sf, ibra_hex_grid) |> as.integer(),
hex_id = ibra_hex_grid$hex_id[intersection]) |>
st_drop_geometry() |>
select(species, hex_id) |>
distinct()
head(species_hex)# A tibble: 6 × 2
species hex_id
<chr> <int>
1 Melithreptus albogularis 836
2 Melithreptus albogularis 1222
3 Melithreptus validirostris 358
4 Melithreptus lunatus 611
5 Melithreptus albogularis 1114
6 Melithreptus chloropsis 629With these unique combinations, we can bring the additional hexagon information into the species_hex data-frame with a left join.
species_hex_combos <- species_hex |>
left_join(ibra_hex_grid, by = "hex_id") |>
st_as_sf()
head(species_hex_combos)Simple feature collection with 6 features and 2 fields
Active geometry column: hex_sf
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 114.7239 ymin: -44.77699 xmax: 153.4064 ymax: -11.77699
Geodetic CRS: GDA94
# A tibble: 6 × 4
species hex_id hex_sf hex_centre_sf
<chr> <int> <POLYGON [°]> <POINT [°]>
1 Melithreptus albog… 836 ((151.097 -26.77699, 151… (152.2517 -26.7766)
2 Melithreptus albog… 1222 ((130.3124 -12.77699, 13… (131.4671 -12.77695)
3 Melithreptus valid… 358 ((145.9009 -43.77699, 14… (147.0556 -43.77521)
4 Melithreptus lunat… 611 ((147.6329 -34.77699, 14… (148.7876 -34.77612)
5 Melithreptus albog… 1114 ((144.1688 -16.77699, 14… (145.3235 -16.7769)
6 Melithreptus chlor… 629 ((114.7239 -33.77699, 11… (115.8786 -33.77619)Species Point Locations
We have now identified the species and hexagon combinations that we wish to visualise. Many hexagons will contain multiple species, and we need a way to show, at minimum, the presence/absence of each species in each hexagon. We will do this using coloured points in a fixed pattern within each hexagon (alternatively you may use differently shaped points).
Given we want to fit a maximum of seven uniquely located points within a hexagon, we can organise them as a smaller hexagon with one point inside each corner, plus a point in the centre [SHOW AN EXAMPLE DESIGN ON THE SIDE/BELOW]. We can then assign each species to one of these seven point locations. It’s worth noting that with a bit of creativity, it is possible to fit any number of points from 2-9 into a hexagon (2-7) or square symmetrically.
You will remember that earlier on, we stored the centre of each hexagon in the data-frame ibra_hex_grid. Each value in the hex_centre_sf column is an sf object of type POINT, as opposed to those in the hex_sf column with are of type POLYGON. With this centre point, we already have one of the seven locations we require for each hexagon. The other six points will be placed halfway between this centre and their vertex.
{sf} has a function that allows us to extract the vertices of each column. st_coordinates() takes the recognised spatial column in an sf data-frame (in our case ibra_hex_grid$hex_sf) and outputs a matrix with a row for each vertex of each feature in the spatial column. The output matrix has four columns - X and Y denote the two coordinates of each point, L1 (which we can ignore), and L2 which contains the row number from ibra_hex_grid that the vertex was taken from.
We can transform this matrix into a tibble and convert the X and Y columns into sf POINTS that each denote a vertex. We then wish to reassign the hex_id column to each vertex, and we can do this by creating a new column with mutate() that indexes the “L2th” hex_id from the initial data-frame. Next, we assign an ID to each type of vertex so that each species point sits in the same relative position across multiple hexagons - since each hexagon is formed in the same way, we can group the data by L2 and assign each an ID based on row-number order. The last step is to remove unneeded columns, rename the new geometry column as hex_vertex_sf and convert the sf data-frame back into a plain tibble. sf dataframes cannot be left/right-joined, but tibbles can, and we will be doing that very soon.
hex_vertices <- ibra_hex_grid |>
st_coordinates() |>
as_tibble() |>
st_as_sf(coords = c("X", "Y"), crs = st_crs(crs), remove = TRUE) |>
group_by(L2) |>
mutate(hex_id = ibra_hex_grid$hex_id[L2],
hex_vertex = row_number()) |>
ungroup() |>
select(-L1, -L2) |>
rename(hex_vertex_sf = geometry) |>
as_tibble()
head(hex_vertices, n = 8)# A tibble: 8 × 3
hex_vertex_sf hex_id hex_vertex
<POINT [°]> <int> <int>
1 (158.0252 -54.77699) 54 1
2 (158.6026 -55.77699) 54 2
3 (159.7573 -55.77699) 54 3
4 (160.3346 -54.77699) 54 4
5 (159.7573 -53.77699) 54 5
6 (158.6026 -53.77699) 54 6
7 (158.0252 -54.77699) 54 7
8 (73.15473 -53.77699) 57 1You may have noticed that each hexagon has seven vertices in hex_vertices, although hexagons only have 6 vertices. Before reading on, think about why this might be. If you’re still unsure, have a look at the coordinates of vertex numbers 1 and 7 for any given hexagon - you’ll notice that they are identical. To create any sf polygon, the string of vertices must end where it began to enclose the shape, so actually we specify all vertices after the first six rows of each hexagon. We will keep vertex 7 though because that row is about to come in handy for our 7th species that will take the point in the middle of the hexagon.
As you can see, we now have two datasets that need to be brought together - species_hex_combos contains the assignments of species to hexagons, as well as sf data on each hexagon; hex_vertices contains the coordinates of each (numbered) vertex of each hexagon.
One task remains before we can join these data-frames - we haven’t assigned our seven species to the seven points in the hexagon. st_make_grid() builds all hexagons in the same way, such that the numbers in our hex_vertex column represent the same vertices, as presented in the diagram below. We are going to provide each species with an ID number that matches the vertex it’ll be represented at (species 7 will go in the middle). We will simply number our species 1 through 7 in the order they appear in species_data, and join this ID column to species_hex_combos. You can choose any order you like.
species_data$species_id <- 1:7
species_hex_combos <- species_hex_combos |>
left_join(species_data |> select(species, species_id), by = "species")
head(species_hex_combos)Simple feature collection with 6 features and 3 fields
Active geometry column: hex_sf
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: 114.7239 ymin: -44.77699 xmax: 153.4064 ymax: -11.77699
Geodetic CRS: GDA94
# A tibble: 6 × 5
species hex_id hex_sf hex_centre_sf species_id
<chr> <int> <POLYGON [°]> <POINT [°]> <int>
1 Melithr… 836 ((151.097 -26.77699, 151… (152.2517 -26.7766) 1
2 Melithr… 1222 ((130.3124 -12.77699, 13… (131.4671 -12.77695) 1
3 Melithr… 358 ((145.9009 -43.77699, 14… (147.0556 -43.77521) 7
4 Melithr… 611 ((147.6329 -34.77699, 14… (148.7876 -34.77612) 2
5 Melithr… 1114 ((144.1688 -16.77699, 14… (145.3235 -16.7769) 1
6 Melithr… 629 ((114.7239 -33.77699, 11… (115.8786 -33.77619) 6We can now join up our two spatial data-frames, and assign species points to locations inside hexagons. We begin by joining hex_vertices to species_hex_combos by matching up rows with the same hex_id and where species_id matches hex_vertex. With spatial columns for both vertex (hex_vertex_sf) and centre (hex_centre_sf), we can create the location we want to place each point point_loc by combining the two columns with st_union() and then calculating their midpoint with st_centroid(). We do this for all species except number 7, which we will just keep the centre instead of calculating a midpoint. All that then remains is to clear the active geometry (st_drop_geometry()), select the important columns for plotting the species points, and reinstate the data-frame as a spatial object (st_as_sf())
species_hex_points <- species_hex_combos |>
left_join(hex_vertices, by = c("hex_id", "species_id" = "hex_vertex")) |>
mutate(point_loc = if_else(species_id == 7,
hex_centre_sf,
st_union(hex_centre_sf, hex_vertex_sf) |> st_centroid())) |>
st_drop_geometry() |>
select(species, hex_id, species_id, point_loc) |>
st_as_sf()All that remains is to plot our data with the set of data-frames we’ve created throughout the process. We need to plot three key features - the background map of Australia, the hexagon outlines, and the coloured species points. We can plot the map of Australia with state boundaries using aus and follow that up with a more detailed outline including island territories with ibra. Distinct hexagon outlines can be obtained by pulling out the unique hexagon geometries from species_hex_combos, and lastly the individual species points can be brought in with species_hex_points and will be coloured by their species argument. Extra colour assignments and customisation can be followed up as well for aesthetic purposes.
hex_map <- ggplot() +
geom_sf(data = aus,
col = "grey20", fill = NA, linewidth = 0.3) +
geom_sf(data = ibra,
col = "grey20", fill = NA, linewidth = 0.1) +
geom_sf(data = species_hex_combos |> select(hex_sf) |> distinct(), aes(geometry = hex_sf),
fill = NA, col = "grey60", linewidth = 0.5) +
geom_sf(data = species_hex_points, aes(geometry = point_loc, col = species),
alpha = 0.95, size = 1.5) +
scale_colour_manual(values = set_names(species_data$colour, species_data$species),
labels = species_data$label) +
lims(x = c(105, 160), y = c(-47, -7)) +
# expand_limits(x = 200) +
guides(colour = guide_legend(
title = element_blank(),
label.theme = element_text(colour = "gray10", size = 8),
byrow = TRUE
)
) +
theme_void() +
theme(legend.justification = c(0, 0),
legend.position = c(0.05, 0.05),
legend.spacing.y = unit(-0.04, "in"),
legend.text = element_markdown())hex_map
Expand for session info
─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.3.2 (2023-10-31)
os macOS Sonoma 14.1.1
system aarch64, darwin20
ui X11
language (EN)
collate en_US.UTF-8
ctype en_US.UTF-8
tz Australia/Melbourne
date 2023-11-27
pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.1)
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.0)
galah * 2.0.0 2023-11-20 [1] CRAN (R 4.3.1)
ggplot2 * 3.4.4 2023-10-12 [1] CRAN (R 4.3.1)
ggtext * 0.1.2 2022-09-16 [1] CRAN (R 4.3.0)
glue * 1.6.2 2022-02-24 [1] CRAN (R 4.3.0)
hexbin * 1.28.3 2023-03-21 [1] CRAN (R 4.3.0)
htmltools * 0.5.7 2023-11-03 [1] CRAN (R 4.3.1)
lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.1)
monochromeR * 0.2.0 2023-09-06 [1] CRAN (R 4.3.0)
ozmaps * 0.4.5 2021-08-03 [1] CRAN (R 4.3.0)
purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.0)
readr * 2.1.4 2023-02-10 [1] CRAN (R 4.3.0)
sessioninfo * 1.2.2 2021-12-06 [1] CRAN (R 4.3.0)
sf * 1.0-14 2023-07-11 [1] CRAN (R 4.3.0)
stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.1)
tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.0)
tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.3.0)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.0)
[1] /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library
──────────────────────────────────────────────────────────────────────────────
